!     Copyright (C) 2010 Benjamin Piaud
!
!     LIMBES is free software; you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation; either version 3, or (at your option)
!     any later version.

!     LIMBES is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.

!     You should have received a copy of the GNU General Public License
!     along with LIMBES; if not, see <http://www.gnu.org/licenses/>

program LIMBES_f_init
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 :: ni,xi,yi,vi
  real*8 :: rho,ux,uy,T,RT,x,y,dx,Lx,Ly
  real*8,allocatable :: feq(:)
  integer*8 :: nv,nvnv
  integer*8 :: n1,n2,n3,n4,n5,n6,n7
  real*8 :: r1,r2,r3,r4,r5,r6,r7,r8
  real*8 :: T0
  character*1 :: char,char1,char2
  character*50 :: arg

  !     pi=3.14
  pi=4.0d+0*datan(one)
  if (iargc().EQ.0) then
     print*,'You forget to indicate the parameter file !!!'
     stop
  endif
  call getarg(1,arg)
  !     fetching simulation parameters
!  open(10,file='../data/simul_parameters.in')
  open(10,file=arg)
  read(10,*) ! time steps number 
  read(10,*) SIMUL_step_number
  read(10,*) ! number of results to save 
  read(10,*) SIMUL_save_number
  read(10,*) ! reference temperature
  read(10,*) T0
  read(10,*) ! gravity
  read(10,*) SIMUL_gx,SIMUL_gy
  read(10,*) ! quadrature order 
  read(10,*) SIMUL_Nquad
  read(10,*) ! Hermite expansion order
  read(10,*) LIMBES_NH
  read(10,*) ! dt/tau ratio
  read(10,*) SIMUL_dttau
  read(10,*) ! transport sheme
  read(10,*) LIMBES_transport_scheme
  read(10,*) ! thermal or isothermal simulation
  read(10,*) SIMUL_thermal
  read(10,*) ! data path
  read(10,*) SIMUL_path
  close(10)

  ! number of discrete velocities
  LIMBES_nv=SIMUL_Nquad**two
  nv=LIMBES_nv
 
  !     fetching mesh information
!  open(10,file='../data/mesh.in')
  open(10,file=SIMUL_path(1:strlen(SIMUL_path))//'mesh.in')
  read(10,*)
  read(10,*) ! cells number
  read(10,*)
  read(10,*) SIMUL_cell_number, SIMUL_dx,SIMUL_nx,SIMUL_ny
  read(10,*) VTK_xmin,VTK_xmax,VTK_ymin,VTK_ymax
  !-------------------------------------------------------------------
  print*,SIMUL_cell_number, SIMUL_dx,SIMUL_nx,SIMUL_ny
  !allocation of LIMBES_var
  call init_LIMBES_var(SIMUL_nx,SIMUL_ny,SIMUL_cell_number,nv,SIMUL_save_number)
  allocate(feq(1:nv))
  !-------------------------------------------------------------------
  do ni=1,SIMUL_cell_number
     read(10,*) n1,n2,n3,n4,n5,n6,n7
     LIMBES_cell_xi(ni)=n2
     LIMBES_cell_yi(ni)=n3
!!$     LIMBES_cell_ni(n2,n3)=n1
!!$     LIMBES_neighbour_RLUD(ni,1)=n4
!!$     LIMBES_neighbour_RLUD(ni,2)=n5
!!$     LIMBES_neighbour_RLUD(ni,3)=n6
!!$     LIMBES_neighbour_RLUD(ni,4)=n7
  enddo
  close(10)

  !     space step evaluation assuming periodic boundary conditions ix x axis
  SIMUL_Lx=SIMUL_dx*dble(SIMUL_nx)
  SIMUL_Ly=SIMUL_dx*dble(SIMUL_ny)
  dx=SIMUL_dx

  !     fetching fluid properties
!  open(10,file='../data/fluid_properties.in')
  open(10,file=SIMUL_path(1:strlen(SIMUL_path))//'fluid_properties.in')
  read(10,*) ! massique fluid constant 
  read(10,*) FLUID_R
  read(10,*) ! collision time
  read(10,*) FLUID_tau
  read(10,*) ! collision time expression
  read(10,*) FLUID_tau_expr
  close(10)

  FLUID_RT0=FLUID_R*T0
  SIMUL_dt=SIMUL_dttau*FLUID_tau
 
  !     initialazing discrete velocities model 
  call LIMBES_dvm()

  open(10,file=SIMUL_path(1:strlen(SIMUL_path))//'f_init.dat',form='unformatted')
 !     intializing macroscopic fields
  do ni=1,SIMUL_cell_number
     xi=LIMBES_cell_xi(ni)
     yi=LIMBES_cell_yi(ni)
     x=dx*dble(xi) - dx/two
     y=dx*dble(yi) - dx/two
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!! macroscopic field !!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     ux=zero            
     uy=zero               
     rho=1.177D+0 - 0.002*x/(5.0d-5)
     T=300.0D+0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     RT=T*FLUID_R
     call LIMBES_calc_feq(rho,ux,uy,RT,feq)
     LIMBES_f(ni,:)=feq(:)
     write(10) (LIMBES_f(ni,vi),vi=1,nv)
  enddo
  close(10)

  SIMUL_dt=SIMUL_dttau*FLUID_tau
  !     few info
  write(*,*) '---------LIMBES init. distribution function---------'
  write(*,*) 'Domain size',SIMUL_nx,'x',SIMUL_ny
  write(*,*) 'CFL=',LIMBES_vx(1)/(SIMUL_dx/SIMUL_dt)

  call free_LIMBES_var()
  deallocate(feq)


end program LIMBES_f_init

